1 SYSTEM SETUP

1.1 Paquetes de R

En primer lugar se cargan los paquetes necesarios para el análisis, y se carga la tabla de datos original “Datos de uso Bicimad Agosto 2018”.

# Load de paquetes para analisis
library(grDevices)
library(geosphere)
library(tidyverse)
library(data.table)
library(plotly)
library(factoextra)
library(ggpubr)
library(Hmisc)
library(corrplot)
library(cluster)
library(NbClust)
library(reprex)
library(magrittr)
library(readr)
library(Rtsne)

# Load dataset target del estudio
# Datos originales de Agosto 2018. 274 mil registros
data_orig <- fread('Bicimad_201808_3.csv', na.strings = '\\N')

2 DATA WRANGLING & FEATURE ENGINEERING

2.1 Tablas auxiliares

Utilizamos datos externos para enriquezer el dataset. Se han utilizado las siguientes fuentes:

Tb. Estaciones - Dataset que contiene el id de estación, dirección, código postal, distrito, coordenadas geográficas… Estos datos se han enriquecido para las estaciones tanto de origen como destino.
Tb. Meteorológica - Datos metereológicos de Madrid por fecha con datos de temperatura y rachas de viento.
Tb. Edad - Tabla que relaciona el código de edad con la descripción literal.
Tb. Tipo de Usuario - Tabla que relaciona el código de tipo de ususario con su descripción literal.

##### LECTURA DE TABLAS #####

# Tb. estaciones -- inicio. Se carga el maestro de estaciones.
master_estaciones_ini <- fread('tb_stations.csv', colClasses = 'character')
master_estaciones_ini <- master_estaciones_ini %>% 
  rename(name_station_ini = name,
         longitude_station_ini = longitude,
         latitude_station_ini = latitude,
         address_station_ini = address,
         idunplug_station_ini = id,
         distrito_station_ini = distrito,
         cp_ini = cod_postal) %>% 
  dplyr::select(-number)


# Tb. estaciones -- fin. 
master_estaciones_fin <- fread('tb_stations.csv', colClasses = 'character')
master_estaciones_fin <- master_estaciones_fin %>% 
  rename(name_station_fin = name,
         longitude_station_fin = longitude,
         latitude_station_fin = latitude,
         address_station_fin = address,
         idplug_station_fin = id,
         distrito_station_fin = distrito,
         cp_fin = cod_postal) %>% 
  dplyr::select(-c(number))

# Tb. Meteo. Se carga y mapea el dataset con los datos climatológicos de agosto.  
master_meteo <-  fread('tb_meteo.csv')
master_meteo <- master_meteo %>%
  rename(unplug_date = fecha) %>%
  mutate(temp_max = as.numeric(sub(",", ".", temp_max, fixed = TRUE)),
         temp_max_hora = as.numeric(sub(",", ".", temp_max_hora, fixed = TRUE)),
         temp_min = as.numeric(sub(",", ".", temp_min, fixed = TRUE)),
         temp_min_hora = as.numeric(sub(",", ".", temp_min_hora, fixed = TRUE)),
         temp_media = as.numeric(sub(",", ".", temp_media, fixed = TRUE)))

# Tb. Edad. Maestro de edad
master_edad <- fread('tb_ageRange.csv', colClasses = 'character')

# Tb. userType. Maestro de tipo de usuarios
master_user <- fread('tb_userType.csv', colClasses = 'character')

En el dataset enriquecido se ha filtrado por código de estación para trabajar únicamente con las estaciones en el distrito “Distrito Retiro”, añadiendo estaciones que por proximidad pueden resultar interesante spara el estudio.

El dataset enriquezido, de ahora en adelante referido como dataset target (var(data)), se ha cruzado por campos comunes y se han normalizados los data types. Además, se ha calculado la distancia mínima entre estaciones computando la Haversine distance entre coordenadas.

# DATA WRANGLING 

# Filtrar por estaciones que nos interesan
estaciones_retiro <- c(64,65,66,69,72,73,74,75,78,79,80,84,85,86,88,90,91,99,100,101,102,107)

# Se unen dataset con los maestros de estaciones, datos meteorológicos, edad y tipo de usuario
data <- data_orig %>% 
  rename(movimiento_id = "id_mov") %>% 
  filter(idunplug_station_ini %in% estaciones_retiro | idplug_station_fin %in% estaciones_retiro) %>% 
  mutate(idunplug_station_ini = as.character(idunplug_station_ini),
         cp_ini = as.character(cp_ini),
         idunplug_base = as.character(idunplug_base),
         idplug_station_fin = as.character(idplug_station_fin),
         cp_fin = as.character(cp_fin),
         idplug_base = as.character(idplug_base),
         user_type = as.character(user_type),
         user_age_range = as.character(user_age_range),
         user_zip_code = as.character(user_zip_code),
         bank_holiday_flg = as.character(bank_holiday_flg),
         unplug_hour = as.character(unplug_hour),
         unplug_date=as.character(unplug_date)
         ) %>% 
  dplyr::select(-c(unplug_hourtime)) %>% 
  left_join(master_edad, by='user_age_range') %>% 
  left_join(master_user, by='user_type') %>% 
  left_join(master_estaciones_ini, by=c('idunplug_station_ini', 'cp_ini')) %>% 
  left_join(master_estaciones_fin, by=c('idplug_station_fin', 'cp_fin')) %>% 
  mutate(longitude_station_ini = as.numeric(longitude_station_ini),
         latitude_station_ini = as.numeric(latitude_station_ini),
         longitude_station_fin = as.numeric(longitude_station_fin),
         latitude_station_fin = as.numeric(latitude_station_fin)) %>% 
  mutate(distance = distHaversine(cbind(longitude_station_ini, latitude_station_ini), cbind(longitude_station_fin, latitude_station_fin))) %>%
           left_join(master_meteo, by='unplug_date')


rm(data_orig, master_edad, master_estaciones_fin, master_estaciones_ini, master_meteo, master_user, estaciones_retiro)

# Ordenar variable 'day'
data$day <- factor(data$day, levels=c('L','M','X','J','V','S','D'))

# save(data, file='data.RData')

3 ANALISIS EXPLORATORIO

El posterior análisis exploratorio describe el dataset target por variables individuales y cruzes entre variables mediante análisis gráficos y sus estadíticos principales.

3.0.1 Distrito de inicio

Más allá de que el distrito 02 - retiro sea el mayoritario (pues estamos influyendo en ello), se observa que los distritos más recurrentes como inicio de viaje son los colindantes al propio retiro, aunque destaca que distrito se situe por encima de Salamanca, cuando es este último el distrito más cercano al propio parque, por lo que ya denota una tendencia natural de los usuarios.

x <- data %>% group_by(distrito_station_ini) %>% summarise(n = length(user_id))
rm(x)

x <- data %>% 
  group_by(distrito_station_ini) %>% 
  summarise(n = length(distrito_station_ini)) %>% 
  ungroup() %>% 
  mutate(porc = n/sum(n)) %>% 
  filter(porc*100 >= 0.1)

z <- ggplot(x, aes(x=distrito_station_ini, y=porc*100))+theme_classic()+
  geom_bar(stat = "identity",fill="lightblue")+ labs(title="Distrito de partida", x="", y="%")+theme(axis.text.x = element_text(angle = 45, hjust = 1))

ggplotly(z)

rm(x, z)

3.0.2 Codigo postal de inicio

Los códigos postales de las estaciones de inicio directamente vinculados con la situación de las estaciones siguen la misma distribución.

x <- data %>% 
  group_by(cp_ini) %>% 
  summarise(n = length(cp_ini)) %>% 
  ungroup() %>% 
  mutate(porc = n/sum(n))

# p <- ggplot(x, aes(x=cp_ini, y=porc*100)) +
#   geom_bar(stat = "identity") +
#   theme(axis.text.x = element_text(angle = 45, hjust = 1))

p <- ggplot(x, aes(x=cp_ini, y=porc*100))+theme_classic()+
  geom_bar(stat = "identity",fill="lightblue")+ labs(title="Código postal de inicio", x="", y="%")+ theme(axis.text.x = element_text(angle = 45, hjust = 1))

ggplotly(p)

rm(x, p)

3.0.3 Estacion de inicio

# 1.2 ID UNPLUG_STATION INI
x <- data %>% group_by(idunplug_station_ini) %>% summarise(n = length(user_id))

x <- data %>% 
  group_by(idunplug_station_ini) %>% 
  summarise(n = length(idunplug_station_ini)) %>% 
  ungroup() %>% 
  mutate(porc = n/sum(n)) %>% 
  filter(porc*100 >= 0.1)

p <- ggplot(x, aes(x=idunplug_station_ini, y=porc*100))+theme_classic()+
   geom_bar(stat = "identity",fill="lightblue")+ labs(title="Estación de partida", x="", y="%") + theme(axis.text.x = element_text(angle = 90))
ggplotly(p)

rm(x,p)

3.0.4 Distrito de fin

Para los distritos en los que finalizan los trayectos, vuelve a aparecer una tendencia en dirección a la zona centro superando a Salamanca que es la distrito más cercano al parque del retiro.

x <- data %>% group_by(distrito_station_fin) %>% summarise(n = length(user_id))
rm(x)

x <- data %>% 
  group_by(distrito_station_fin) %>% 
  summarise(n = length(distrito_station_fin)) %>% 
  ungroup() %>% 
  mutate(porc = n/sum(n)) %>% 
  filter(porc*100 >= 0.1)

z <- ggplot(x, aes(x=distrito_station_fin, y=porc*100))+theme_classic()+
  geom_bar(stat = "identity",fill="lightblue")+ labs(title="Distrito de llegada", x="", y="%")+ theme(axis.text.x = element_text(angle = 45, hjust = 1))

ggplotly(z)

rm(x, z)

3.0.5 Codigo postal de fin

Distribución calcada con las estaciones de inicio, pues estamos eligiendo aquellos viajes que comienzan o terminan en el retiro.

x <- data %>% 
  group_by(cp_fin) %>% 
  summarise(n = length(cp_fin)) %>% 
  ungroup() %>% 
  mutate(porc = n/sum(n))

# p <- ggplot(x, aes(x=cp_fin, y=porc*100)) + 
#   geom_bar(stat = "identity") + 
#   theme(axis.text.x = element_text(angle = 45, hjust = 1))

p <- ggplot(x, aes(x=cp_fin, y=porc*100))+theme_classic()+
    geom_bar(stat = "identity",fill="lightblue")+ labs(title="Codigo postal de fin", x="", y="%")
  ggplotly(p)

rm(x, p)

3.0.6 Estacion de fin

Conclusiones similares a las estaciones de inicio, donde al estar filtrando por los viajes que se inician o terminan en el retiro estamos forzando que estas estaciones sean más frecuentes.

x <- data %>% 
  group_by(idplug_station_fin) %>% 
  summarise(n = length(idplug_station_fin)) %>% 
  ungroup() %>% 
  mutate(porc = n/sum(n)) %>% 
  filter(porc*100 >= 0.1)

# p <- ggplot(x, aes(x=idplug_station_fin, y=porc*100)) +
#   geom_bar(stat = "identity") +
#   theme(axis.text.x = element_text(angle = 45, hjust = 1))

p <- ggplot(x, aes(x=idplug_station_fin, y=porc*100))+theme_classic()+
  geom_bar(stat = "identity",fill="lightblue")+ labs(title="Estación de llegada", x="", y="%")

ggplotly(p)

rm(x, p)

3.0.7 Tipo de usuario

Los usuarios conabono anual son muy claramente mayoritarios, estando presentes en el 86% de los viajes que se inician o terminan en el retiro, con menos de un 10% los trabajadores de Bicimad y con el 4% los usuarios con abono ocasional.

x <- data %>% 
  group_by(user_type_lit) %>% 
  summarise(n = length(user_type_lit)) %>% 
  ungroup() %>% 
  mutate(porc = n/sum(n))
  

# p <- ggplot(x, aes(x=user_type, y=porc*100)) + 
#   geom_bar(stat = "identity") + 
#   theme(axis.text.x = element_text(angle = 45, hjust = 1))

p <- ggplot(x, aes(x=user_type_lit, y=porc*100))+theme_classic()+
    geom_bar(stat = "identity",fill="lightblue")+ labs(title="Tipo de usuario", x="", y="%")

ggplotly(p)

rm(x, p)

3.0.8 Rango de edad de usuario

Anque la mayoría de usuarios no declaran la edad, de aquellos que la declaran el rango. Si excluimos a estos y recalculamos los procentajes estaríamos diciendo que la mitad de los viajes con origen o destino en el retiro son por usuarios de entre 27 y 40 años, donde si unimo los usuarios de 27 a 40, tendríamos un peso del 86% de usuarios adultos. Con esto coge peso la teoría de que es un servicio usado por gente adulta, pudiendo ser usada para desplazamientos en dirección al trabajo, por ejemplo.

x <- data %>% 
  group_by(user_age_range_lit) %>% 
  summarise(n = length(user_age_range_lit)) %>% 
  ungroup() %>% 
  mutate(porc = n/sum(n)) 
  
p <- ggplot(x, aes(x=user_age_range_lit, y=porc*100))  + 
    theme_classic() +
    geom_bar(stat = "identity",fill="lightblue") + 
    labs(title="Rango de edad de usuario", x="Grupo", y="%")+
    theme(axis.text.x = element_text(angle = 45, hjust = 1))
ggplotly(p)

rm(x, p)

3.0.9 Codigo postal de usuario

Aunque gran parte de los usuarios no declaran su codigo portal, en el histograma destaca la perdida de peso del código postal 28014, que si estaba muy presente en las estaciones de inicio y fin. Esto podría ser acausa de que esta localización es más turistica que residencial.

x <-  filter(data,!is.na(user_zip_code))
x <- x %>% 
  group_by(user_zip_code) %>% 
  summarise(n = length(user_zip_code)) %>% 
  ungroup() %>% 
  mutate(porc = n/sum(n))%>% 
  filter(porc*100 >= 0.1)

p <- ggplot(x, aes(x=user_zip_code, y=porc*100))  +
    theme_classic() +
    geom_bar(stat = "identity",fill="lightblue") +
    labs(title="Codigo postal de usuario", x="", y="%")+
   theme(axis.text.x = element_text(angle = 90, hjust = 1))

ggplotly(p)

rm(x, p)

3.0.10 Dia de la semana del trayecto

En base al histograma, el servicio de bicimad para los viajes relativos a retiro, parece ser usado como media de transporte rutinaria, más que a modo ocasional, con una tendencia muy marcada durante la semana, que cae en el fin de semana. Es posible que si cogiesemos un dataset de otro mes que no sea de un periodo vacacional, estas diferencias sean más marcadas todabía.

x <- data %>% 
  group_by(day) %>% 
  summarise(n = length(day)) %>% 
  ungroup() %>% 
  mutate(porc = n/sum(n))


p <- ggplot(x, aes(x=day, y=porc*100))  +
    theme_classic() +
    geom_bar(stat = "identity",fill="lightblue") +
    labs(title="Día de la semana", x="", y="%")

ggplotly(p)

rm(x, p)

3.0.11 Festivos/Laborables

Se observa que el 77% de los trayectos ocurren en dias laborables, y el restate en festivos(incluyendo el dia 15 de Agosto).

x <- data %>% 
  group_by(bank_holiday_flg) %>% 
  summarise(n = length(bank_holiday_flg)) %>% 
  ungroup() %>% 
  mutate(porc = n/sum(n))

p <- ggplot(x, aes(x=bank_holiday_flg, y=porc*100)) + 
     theme_classic() +
     geom_bar(stat = "identity",fill="lightblue") +
     labs(title=" Laborables vs Festivos", x="", y="%")
ggplotly(p)

rm(x, p)

3.0.12 Hora de inicio del trayecto

Se notan horas pico coincidiendo con el inicio (8) y fin (19) de la que podría ser una jornada laboral, acentuando además las 15h entendiendolo como la hora de comer. En un análisis más profundo, filtrando por bank_holiday_flg = 0, es decir, dias laborables, esta diferencia se ve acentuada.

x <- data %>% 
  group_by(unplug_hour,bank_holiday_flg) %>% 
  summarise(n = length(unplug_hour)) %>% 
  ungroup() %>% 
  mutate(porc = n/sum(n))
x$unplug_hour=as.numeric(x$unplug_hour)
p <- ggplot(x, aes(x=unplug_hour, y=porc*100)) + 
    theme_classic() +
    geom_bar(stat = "identity",fill="lightblue") +
    labs(title="Hora de inicio del trayecto", x="", y="%") +
    theme(axis.text.x = element_text(angle = 45, hjust = 1))
ggplotly(p)

rm(x, p)

3.1 ESTADISTICOS PRINCIPALES

3.1.1 Duracion del trayecto

Transformando la variale numérica de travel_time a minutos para hacerla más legible, observamos que aunque existen outliers, con viajes de 8.000 min o de 0.17 min, la media de los viajes está entorno a 11 minutos, situándose el 3er cuartil por debajo de los 20 minutos.

data <- data %>% 
  mutate(travel_time_mins = travel_time/60)

Descripcion de estadisticos:

summary(data$travel_time_mins)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##    0.017    7.817   11.400   20.827   18.050 8030.233

Desviacion estandar:

sd(data$travel_time_mins)

## [1] 87.87702

Varianza:

var(data$travel_time_mins)

## [1] 7722.371

3.1.2 Distancia del trayecto

Observamos que la distancia media recorrida está entorno a los 2KM. La mayoría de los viajes están contenidos en la franja de los 1-3km.

Descripcion de estadisticos:

summary(data$distance)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0    1235    1823    1889    2397    5944

Desviacion tipica:

sd(data$distance)

## [1] 956.4534

Varianza:

var(data$distance)

## [1] 914803.1

3.1.3 Temperatura maxima

Altas temperaturas sin bajas de los 30 grados, sin haber grandes dispersiones. De haber sido otros meses para analizar, habría cogido valor, el introducir datos referentes a las precipitaciones

Descripcion de estadisticos:

summary(data$temp_max)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   30.40   32.70   33.80   34.43   36.70   40.00

Desviacion tipica:

sd(data$temp_max)

## [1] 2.508407

Varianza:

var(data$temp_max)

## [1] 6.292106

3.1.4 Temperatura minima

De igual manera las temperaturas minimas con poca dispersión e igualmente altas.

Descripcion de estadisticos:

summary(data$temp_min)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   16.80   19.60   21.00   21.36   23.00   25.90

Desviacion tipica:

sd(data$temp_min)

## [1] 2.341467

Varianza:

var(data$temp_min)

## [1] 5.48247

3.1.5 Velocidad del viento maximas

Descripcion de estadisticos:

summary(data$viento_max)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   11.00   13.00   15.00   16.26   19.00   25.00

Desviacion tipica:

sd(data$viento_max)

## [1] 3.597168

Varianza:

var(data$viento_max)

## [1] 12.93962

4 ANALISIS DE OUTLIERS

4.1 TRAVEL_TIME

En primera instancia, cuando tratamos de analizar la variable travel time, los outliers no nos dejan ver realmente como se comportan los usuarios, en el momento que los dejamos de lado para pintar el boxplot, observamos lo que nos desprendió el análisis estadistico, donde el servicio era usado para viajes relativamente cortos en su mayoría.

# 1.C.1 TRAVEL_TIME
# Anadir 'travel_time' en minutos
data <- data %>% 
  mutate(travel_time_mins = travel_time/60)

# Identificacion de outliers
outliers_travel_time_mins <- boxplot.stats(data$travel_time_mins)$out
p <- ggplot(data, aes(y=travel_time_mins)) + geom_boxplot() + labs(title="Travel_time - outliers")
ggplotly(p)

# Quitando outliers de 'data'
data_no_outliers <- data[-which(data$travel_time_mins %in% outliers_travel_time_mins),]
# save(data_no_outliers, file='data_no_outliers.RData')

# Boxplot - 'data_no_outliers'
p <- ggplot(data_no_outliers, aes(y=travel_time_mins)) + geom_boxplot() + labs(title="Travel_time - no outliers")
ggplotly(p)

rm(p, outliers_travel_time_mins)

4.2 DISTANCE

Los outliers no dificultan tanto el análisis de la variable en este caso, aun así en el momento que los excluimos observamos que en terminos generales los viajes rondan los 1 a 3 KM.

# 1.C.2 DISTANCE
# Identificacion de outliers
outliers_distance <- boxplot.stats(data$distance)$out
p <- ggplot(data, aes(y=distance)) + geom_boxplot() + labs(title="Distance - outliers")
ggplotly(p)

# Quitando outliers de 'data'
data_no_outliers_mod <- data[-which(data$distance %in% outliers_distance),]

# Boxplot - 'data_no_outliers'
p <- ggplot(data_no_outliers_mod, aes(y=distance)) + geom_boxplot() + labs(title="Distance - no outliers")
ggplotly(p)

rm(p, outliers_distance, data_no_outliers_mod)

5 ANALISIS BIVARIANTE

Para conocer la relación entre dos variables categóricas utilizaremos:

Tablas de contingencia para describir esta relación a nivel muestral. Test de chi-cuadrado para hacer inferencia sobre la relación entre ambas variables.

5.1 Tb.Contingencia - Edad de usuario v. Tipo de usuario

El test de Chi cuadrado no muestra que haya una correlación clara.
Hay 2 rangos principales (27-65 años) que aglutinan la mayor parte de valores. Esto nos llevará a hacer un análisis más profundo para determinar si el uso principal de bicimad son trabajadores.
Hay muchos usuarios con edad desconocida lo que desvirtúa el dataset.
Los usuarios ocasionales no suelen registrar su edad, por lo que podemos inferir que su representatividad porcentual dentro del grupo “unknown” es más elevada que para trabajadores o usuarios anuales.

#0 - Unknown
#1 - 0 y 16 años
#2 - 17 y 18 años
#3 - 19 y 26 años
#4 - 27 y 40 años
#5 - 41 y 65 años
#6 - > 66 años

frec <- round(prop.table(table(data$user_type, data$user_age_range))*100, 2)
frec

##    
##         0     1     2     3     4     5     6
##   1 34.83  0.20  0.71  5.56 27.53 17.37  0.50
##   2  3.83  0.00  0.00  0.02  0.01  0.01  0.00
##   3  3.49  0.74  0.00  0.28  2.34  2.59  0.00

#Calculamos "Chi-cuadrado"
chisq.test(table(data$user_type, data$user_age_range))

## 
##  Pearson's Chi-squared test
## 
## data:  table(data$user_type, data$user_age_range)
## X-squared = 6806.9, df = 12, p-value < 2.2e-16

#Visualizamos la relación en un histograma
#heatmap(table(data$user_type, data$user_age_range), Colv = NA, Rowv = NA, xlab="user_age_range", ylab="user_type")

barplot(frec, beside=FALSE,legend.text= rownames(frec[1]), main='Distribución % Tipo usuario por edad', xlab='Rangos Edad', ylab='% Distribución tipo usuario')

rm(frec)

5.2 Estudio por Edad de los trayectos con distrito final enel Retiro y distrito inicial diferente

Por el valor de P-value se concluye que no haya correlación entre variables.
Destacan el distrito Centro y Salamanca como distritos con más viajes de inicio
Tetuan y Moncloa los que menos viajes presentan.
Respecto a los rangos de Edad tiene lógica que los Jubilados que viven en Salamanca sean los que más vayan al retiro.

#Generamos la tabla de variables eliminando como distrito inicial el Retiro 
table_distrit_age <- data %>% 
  dplyr::select(c(distrito_station_ini,user_age_range)) %>% 
  filter(grepl("03",distrito_station_ini)==FALSE) %>% 
  mutate(distrito_station_ini = as.factor(distrito_station_ini),
         user_age_range = as.factor(user_age_range))

#Realizamos la matriz de contingencia de las frecuencias relativas
round(prop.table(table(table_distrit_age))*100, 2)

##                     user_age_range
## distrito_station_ini     0     1     2     3     4     5     6
##       01  CENTRO     17.04  0.21  0.25  2.17 11.90  6.68  0.10
##       02  ARGANZUELA  7.56  0.25  0.13  1.06  4.93  3.28  0.04
##       04  SALAMANCA  13.01  0.26  0.19  1.64  9.57  6.05  0.21
##       05  CHAMARTÍN   1.01  0.04  0.05  0.24  0.88  0.69  0.00
##       06  TETUÁN      0.52  0.02  0.00  0.06  0.39  0.24  0.00
##       07  CHAMBERÍ    3.00  0.04  0.04  0.47  2.11  1.59  0.03
##       09  MONCLOA     0.86  0.00  0.02  0.08  0.65  0.43  0.01

#Incluimos simulación Chi-cuadrado 
chisq.test(table_distrit_age$distrito_station_ini,table_distrit_age$user_age_range)

## 
##  Pearson's Chi-squared test
## 
## data:  table_distrit_age$distrito_station_ini and table_distrit_age$user_age_range
## X-squared = 215.74, df = 36, p-value < 2.2e-16

#Obtenemos un Warning al tener valores inferiores a 5 y realizamos también el Test de fisher
fisher.test(table_distrit_age$distrito_station_ini,table_distrit_age$user_age_range,simulate.p.value=TRUE)

## 
##  Fisher's Exact Test for Count Data with simulated p-value (based
##  on 2000 replicates)
## 
## data:  table_distrit_age$distrito_station_ini and table_distrit_age$user_age_range
## p-value = 0.0004998
## alternative hypothesis: two.sided

#Visualizamos la relación en un mapa de calor
heatmap( xtabs(~  user_age_range + distrito_station_ini, table_distrit_age), ylab="user_age_range",margins =c(10,7))

rm(table_distrit_age)

5.3 Estudio por Edad de los trayectos con origen distrito Retiro y distrito final diferente al Retiro.

El test de Fisher devuelve un concluye que no hay correlación entre variables.
Destacan el distrito Centro y Retiro como distritos con mayor número de trayectos como destino final.
Moncloa y Chamartin son los que menos movimientos presentan cómo destino final para trayectos iniciados en el Retiro.

#Extraemos la variable distrito inicial con movimientos que parte solo del Retiro 
table_distrit_fin_age<- data %>% 
  dplyr::select(c(distrito_station_fin,user_age_range,distrito_station_ini)) %>% 
  filter(grepl("03",distrito_station_ini)==TRUE)%>%
  mutate(distrito_station_fin=as.factor(distrito_station_fin),
         user_age_range=as.factor(user_age_range))%>%
  dplyr::select(-c(distrito_station_ini))

#Matriz de contingencia destacando frecuencia de edad (Sumatoria cols = 1)
round(prop.table(table(table_distrit_fin_age))*100, 2)

##                     user_age_range
## distrito_station_fin     0     1     2     3     4     5     6
##       01  CENTRO     11.87  0.11  0.15  1.94  8.64  5.29  0.09
##       02  ARGANZUELA  3.66  0.09  0.07  0.44  2.32  1.76  0.01
##       03  RETIRO     14.04  0.39  0.31  2.15  8.54  8.28  0.37
##       04  SALAMANCA   6.97  0.26  0.11  0.90  5.88  3.45  0.16
##       05  CHAMARTÍN   0.95  0.05  0.01  0.15  0.87  0.58  0.00
##       06  TETUÁN      0.45  0.08  0.00  0.10  0.22  0.22  0.01
##       07  CHAMBERÍ    2.38  0.10  0.06  0.26  2.13  1.20  0.02
##       09  MONCLOA     0.79  0.00  0.03  0.08  0.59  0.45  0.00

#Realizamos también el Test de fisher
fisher.test(table(table_distrit_fin_age),simulate.p.value=TRUE)

## 
##  Fisher's Exact Test for Count Data with simulated p-value (based
##  on 2000 replicates)
## 
## data:  table(table_distrit_fin_age)
## p-value = 0.0004998
## alternative hypothesis: two.sided

#Visualizamos la relación 
heatmap( xtabs(~  user_age_range + distrito_station_fin, table_distrit_fin_age), ylab="user_age_range",margins =c(10,7))

rm(table_distrit_fin_age)

6 ANALISIS DE CORRELACION

El objetivo es detectar posibles relaciones lineales entre las variables. Si presentan un fuerte grado de correlación, buscar la forma funcional que mejor explique la variable dependiente a partir de la independiente (Análisis de regresión)

# Dataset con solo variables continuas
data_cont <- select_if(data, is.numeric)
summary(data_cont)

##   travel_time     longitude_station_ini latitude_station_ini
##  Min.   :     1   Min.   :-3.725        Min.   :40.39       
##  1st Qu.:   469   1st Qu.:-3.694        1st Qu.:40.41       
##  Median :   684   Median :-3.688        Median :40.42       
##  Mean   :  1250   Mean   :-3.688        Mean   :40.42       
##  3rd Qu.:  1083   3rd Qu.:-3.678        3rd Qu.:40.42       
##  Max.   :481814   Max.   :-3.669        Max.   :40.46       
##  longitude_station_fin latitude_station_fin    distance       temp_max    
##  Min.   :-3.725        Min.   :40.39        Min.   :   0   Min.   :30.40  
##  1st Qu.:-3.694        1st Qu.:40.41        1st Qu.:1235   1st Qu.:32.70  
##  Median :-3.688        Median :40.42        Median :1823   Median :33.80  
##  Mean   :-3.688        Mean   :40.42        Mean   :1889   Mean   :34.43  
##  3rd Qu.:-3.678        3rd Qu.:40.42        3rd Qu.:2397   3rd Qu.:36.70  
##  Max.   :-3.669        Max.   :40.46        Max.   :5944   Max.   :40.00  
##  temp_max_hora      temp_min     temp_min_hora     temp_media   
##  Min.   :15.00   Min.   :16.80   Min.   :5.000   Min.   :23.60  
##  1st Qu.:15.00   1st Qu.:19.60   1st Qu.:7.000   1st Qu.:26.10  
##  Median :16.00   Median :21.00   Median :8.000   Median :27.50  
##  Mean   :15.84   Mean   :21.36   Mean   :7.316   Mean   :27.89  
##  3rd Qu.:16.00   3rd Qu.:23.00   3rd Qu.:8.000   3rd Qu.:30.10  
##  Max.   :18.00   Max.   :25.90   Max.   :8.000   Max.   :32.90  
##    viento_max    viento_max_hora travel_time_mins  
##  Min.   :11.00   Min.   : 0.00   Min.   :   0.017  
##  1st Qu.:13.00   1st Qu.: 3.00   1st Qu.:   7.817  
##  Median :15.00   Median :18.00   Median :  11.400  
##  Mean   :16.26   Mean   :13.26   Mean   :  20.827  
##  3rd Qu.:19.00   3rd Qu.:19.00   3rd Qu.:  18.050  
##  Max.   :25.00   Max.   :23.00   Max.   :8030.233

Generamos la matriz de coeficientes de correlación de Pearson y Spearman para todos los pares de variables:

Correlación de Pearson (r): tendencias lineales
Correlación de Spearman (s): tendencias monótonas (ritmo no constante)

6.1 Matriz de correlación - Pearson

La geolocalización de las estaciones guarda ligera relación con la distancia.
Las variables meteorológicas tienen una leve correlación.

# Matriz de correlacion de Pearson
corr_matrix_p <- rcorr(as.matrix(data_cont))

# Mostramos los valores de los coeficientes y no se identifica ningún valor cercano a 1 o -1
corr_matrix_p[["r"]]

##                         travel_time longitude_station_ini
## travel_time            1.0000000000          -0.017674873
## longitude_station_ini -0.0176748728           1.000000000
## latitude_station_ini  -0.0046979295          -0.019512680
## longitude_station_fin -0.0206480997          -0.097530292
## latitude_station_fin   0.0127003206           0.060140530
## distance               0.0274774723          -0.151059373
## temp_max              -0.0066242061           0.006488413
## temp_max_hora          0.0093248581          -0.002565709
## temp_min              -0.0064799604           0.015211328
## temp_min_hora         -0.0021280376          -0.002139069
## temp_media            -0.0068135284           0.011269486
## viento_max            -0.0007073532          -0.002400726
## viento_max_hora       -0.0019142138           0.013846186
## travel_time_mins       1.0000000000          -0.017674873
##                       latitude_station_ini longitude_station_fin
## travel_time                   -0.004697929         -0.0206480997
## longitude_station_ini         -0.019512680         -0.0975302917
## latitude_station_ini           1.000000000          0.0577037050
## longitude_station_fin          0.057703705          1.0000000000
## latitude_station_fin           0.004058602         -0.0372414132
## distance                       0.233569721         -0.1508242916
## temp_max                       0.003365522         -0.0003444813
## temp_max_hora                  0.004495453          0.0021605950
## temp_min                       0.008056428          0.0043823536
## temp_min_hora                  0.001616876         -0.0023939381
## temp_media                     0.005966527          0.0019966595
## viento_max                    -0.002656509         -0.0042835918
## viento_max_hora                0.006068944          0.0094009747
## travel_time_mins              -0.004697929         -0.0206480997
##                       latitude_station_fin      distance      temp_max
## travel_time                   0.0127003206  0.0274774723 -0.0066242061
## longitude_station_ini         0.0601405295 -0.1510593735  0.0064884129
## latitude_station_ini          0.0040586023  0.2335697210  0.0033655215
## longitude_station_fin        -0.0372414132 -0.1508242916 -0.0003444813
## latitude_station_fin          1.0000000000  0.2506925000  0.0002142818
## distance                      0.2506925000  1.0000000000 -0.0065910863
## temp_max                      0.0002142818 -0.0065910863  1.0000000000
## temp_max_hora                 0.0051896320 -0.0100235447  0.2798889888
## temp_min                      0.0063245961  0.0004258983  0.7861133440
## temp_min_hora                -0.0095104846  0.0004367628 -0.1289228321
## temp_media                    0.0034923128 -0.0033068938  0.9481690531
## viento_max                   -0.0141966883  0.0038064569 -0.2375628812
## viento_max_hora               0.0096913922  0.0176732898 -0.0548403037
## travel_time_mins              0.0127003206  0.0274774723 -0.0066242061
##                       temp_max_hora      temp_min temp_min_hora
## travel_time             0.009324858 -0.0064799604 -0.0021280376
## longitude_station_ini  -0.002565709  0.0152113281 -0.0021390692
## latitude_station_ini    0.004495453  0.0080564284  0.0016168764
## longitude_station_fin   0.002160595  0.0043823536 -0.0023939381
## latitude_station_fin    0.005189632  0.0063245961 -0.0095104846
## distance               -0.010023545  0.0004258983  0.0004367628
## temp_max                0.279888989  0.7861133440 -0.1289228321
## temp_max_hora           1.000000000  0.0750279766  0.0305656189
## temp_min                0.075027977  1.0000000000 -0.1693298519
## temp_min_hora           0.030565619 -0.1693298519  1.0000000000
## temp_media              0.191772160  0.9416087152 -0.1555022294
## viento_max              0.026212462 -0.2475639319 -0.2612309556
## viento_max_hora        -0.029592376  0.2267748966  0.0256348279
## travel_time_mins        0.009324858 -0.0064799604 -0.0021280376
##                         temp_media    viento_max viento_max_hora
## travel_time           -0.006813528 -0.0007073532    -0.001914214
## longitude_station_ini  0.011269486 -0.0024007262     0.013846186
## latitude_station_ini   0.005966527 -0.0026565093     0.006068944
## longitude_station_fin  0.001996659 -0.0042835918     0.009400975
## latitude_station_fin   0.003492313 -0.0141966883     0.009691392
## distance              -0.003306894  0.0038064569     0.017673290
## temp_max               0.948169053 -0.2375628812    -0.054840304
## temp_max_hora          0.191772160  0.0262124623    -0.029592376
## temp_min               0.941608715 -0.2475639319     0.226774897
## temp_min_hora         -0.155502229 -0.2612309556     0.025634828
## temp_media             1.000000000 -0.2556799958     0.089117332
## viento_max            -0.255679996  1.0000000000     0.040719742
## viento_max_hora        0.089117332  0.0407197419     1.000000000
## travel_time_mins      -0.006813528 -0.0007073532    -0.001914214
##                       travel_time_mins
## travel_time               1.0000000000
## longitude_station_ini    -0.0176748728
## latitude_station_ini     -0.0046979295
## longitude_station_fin    -0.0206480997
## latitude_station_fin      0.0127003206
## distance                  0.0274774723
## temp_max                 -0.0066242061
## temp_max_hora             0.0093248581
## temp_min                 -0.0064799604
## temp_min_hora            -0.0021280376
## temp_media               -0.0068135284
## viento_max               -0.0007073532
## viento_max_hora          -0.0019142138
## travel_time_mins          1.0000000000

# Heatmap Pearson
corrplot(corr_matrix_p$r, type = "upper", order = "alphabet", tl.col = "black", win.asp = 0.8, tl.cex = 0.8 , title="Heatmap - Pearson", mar=c(0,0,1,0))

rm(corr_matrix_p)

6.2 Matriz de correlación - Spearman

Se puede apreciar una ligera correlación directa entre tiempo y distancia

# Matriz de correlación de Spearman
corr_matrix_s <- rcorr(as.matrix(data_cont), type="spearman")

# Mostramos los valores de los coeficientes y no se identifica ningún valor cercano a 1 o -1
corr_matrix_s[["r"]]

##                        travel_time longitude_station_ini
## travel_time            1.000000000          -0.094961698
## longitude_station_ini -0.094961698           1.000000000
## latitude_station_ini   0.035932870           0.024231486
## longitude_station_fin -0.102620135          -0.073095437
## latitude_station_fin   0.106031785           0.040563783
## distance               0.469928266          -0.117000020
## temp_max              -0.032810472           0.003177315
## temp_max_hora          0.006815319          -0.001746893
## temp_min              -0.040524887           0.011266963
## temp_min_hora         -0.006415498           0.001406170
## temp_media            -0.039473714           0.008964384
## viento_max             0.011486869          -0.001551442
## viento_max_hora       -0.025575981           0.009992054
## travel_time_mins       1.000000000          -0.094961698
##                       latitude_station_ini longitude_station_fin
## travel_time                    0.035932870         -0.1026201352
## longitude_station_ini          0.024231486         -0.0730954373
## latitude_station_ini           1.000000000          0.0381286270
## longitude_station_fin          0.038128627          1.0000000000
## latitude_station_fin           0.018471512          0.0039950635
## distance                       0.095878553         -0.1183316867
## temp_max                       0.006916440         -0.0004112203
## temp_max_hora                  0.004927458          0.0026335360
## temp_min                       0.010850878          0.0020978586
## temp_min_hora                  0.005432258         -0.0001942805
## temp_media                     0.010821959          0.0020264022
## viento_max                    -0.006910896         -0.0020979001
## viento_max_hora                0.004396632          0.0078291242
## travel_time_mins               0.035932870         -0.1026201352
##                       latitude_station_fin      distance      temp_max
## travel_time                    0.106031785  0.4699282661 -0.0328104722
## longitude_station_ini          0.040563783 -0.1170000197  0.0031773152
## latitude_station_ini           0.018471512  0.0958785529  0.0069164404
## longitude_station_fin          0.003995064 -0.1183316867 -0.0004112203
## latitude_station_fin           1.000000000  0.1145705635  0.0079022528
## distance                       0.114570564  1.0000000000 -0.0032922756
## temp_max                       0.007902253 -0.0032922756  1.0000000000
## temp_max_hora                  0.007735662 -0.0122418514  0.1528658527
## temp_min                       0.009723055 -0.0010792157  0.7655483162
## temp_min_hora                 -0.003054865  0.0008732965 -0.0461419882
## temp_media                     0.009662459 -0.0038717199  0.9303800647
## viento_max                    -0.009691485  0.0040489509 -0.2393712274
## viento_max_hora                0.011372480  0.0137994236 -0.0805360277
## travel_time_mins               0.106031785  0.4699282661 -0.0328104722
##                       temp_max_hora     temp_min temp_min_hora
## travel_time             0.006815319 -0.040524887 -0.0064154981
## longitude_station_ini  -0.001746893  0.011266963  0.0014061699
## latitude_station_ini    0.004927458  0.010850878  0.0054322583
## longitude_station_fin   0.002633536  0.002097859 -0.0001942805
## latitude_station_fin    0.007735662  0.009723055 -0.0030548652
## distance               -0.012241851 -0.001079216  0.0008732965
## temp_max                0.152865853  0.765548316 -0.0461419882
## temp_max_hora           1.000000000  0.106337542  0.1566595557
## temp_min                0.106337542  1.000000000  0.0112661711
## temp_min_hora           0.156659556  0.011266171  1.0000000000
## temp_media              0.139980278  0.931195119 -0.0059619988
## viento_max              0.025417408 -0.218518385 -0.2680092973
## viento_max_hora         0.111698285  0.088291250  0.0274262549
## travel_time_mins        0.006815319 -0.040524887 -0.0064154981
##                         temp_media   viento_max viento_max_hora
## travel_time           -0.039473714  0.011486869    -0.025575981
## longitude_station_ini  0.008964384 -0.001551442     0.009992054
## latitude_station_ini   0.010821959 -0.006910896     0.004396632
## longitude_station_fin  0.002026402 -0.002097900     0.007829124
## latitude_station_fin   0.009662459 -0.009691485     0.011372480
## distance              -0.003871720  0.004048951     0.013799424
## temp_max               0.930380065 -0.239371227    -0.080536028
## temp_max_hora          0.139980278  0.025417408     0.111698285
## temp_min               0.931195119 -0.218518385     0.088291250
## temp_min_hora         -0.005961999 -0.268009297     0.027426255
## temp_media             1.000000000 -0.226439217    -0.012671728
## viento_max            -0.226439217  1.000000000     0.176807966
## viento_max_hora       -0.012671728  0.176807966     1.000000000
## travel_time_mins      -0.039473714  0.011486869    -0.025575981
##                       travel_time_mins
## travel_time                1.000000000
## longitude_station_ini     -0.094961698
## latitude_station_ini       0.035932870
## longitude_station_fin     -0.102620135
## latitude_station_fin       0.106031785
## distance                   0.469928266
## temp_max                  -0.032810472
## temp_max_hora              0.006815319
## temp_min                  -0.040524887
## temp_min_hora             -0.006415498
## temp_media                -0.039473714
## viento_max                 0.011486869
## viento_max_hora           -0.025575981
## travel_time_mins           1.000000000

# Heatmap Spearman
corrplot(corr_matrix_s$r, type = "upper", order = "alphabet", tl.col = "black", win.asp = 0.8, tl.cex = 0.8, title="Heatmap - Spearman", mar=c(0,0,1,0))

rm(corr_matrix_s)

6.3 Análisis gráfico - distance vs travel_time_mins

ggplot( data_cont , aes(x = distance, y = travel_time_mins)) +
  geom_point() + 
  theme(plot.subtitle = element_text(vjust = 1), 
        plot.caption = element_text(vjust = 1),
        panel.grid.major = element_line(colour = "orange"), 
        axis.title = element_text(face = "bold"), 
        axis.text = element_text(face = "bold"), 
        plot.title = element_text(face = "bold", colour = "chocolate")) +
  labs(title = "Estudio Correlación Distancia - Tiempo Viaje", subtitle = "")

Eliminamos Outliers que como vimos antes son duraciones superiores a 1,30 horas segundos y descartando a los trabajadores

y <- data %>% 
  select(distance, travel_time_mins,user_type) %>% 
  filter(travel_time_mins < 87,
         distance > 0,
         user_type=='1'|user_type=='2') %>%
  dplyr::select(-c(user_type))

ggplot( y , aes(x = distance, y = travel_time_mins)) + 
  geom_point() + 
  geom_smooth(method=lm) + 
  theme(plot.subtitle = element_text(vjust = 1), 
        plot.caption = element_text(vjust = 1), 
        panel.grid.major = element_line(colour = "orange"), 
        axis.title = element_text(face = "bold"), 
        axis.text = element_text(face = "bold"), 
        plot.title = element_text(face = "bold", colour = "chocolate")) + 
  labs(title = "Estudio Correlación Distancia - Tiempo Viaje", subtitle = "")

rm(y)

Se puede apreciar una ligerísima correlación directa entre tiempo y distancia
La hora de viento máximo tiene una muy leve correlación inversa con la temperatura mínima
La geolocalización de las estaciones guarda ligera relación con la distancia.

6.4 Análisis de correlación - user_type vs travel_time_mins

#Extraemos la variable distrito inicial eliminado Retiro 
table_travel_user <- data %>% 
  dplyr::select(c(travel_time_mins, user_type)) %>% 
  mutate(user_type = as.numeric(user_type))

#calculamos coeficiente de pearson
cor.test(table_travel_user$travel_time_mins, table_travel_user$user_type, method = c("pearson"))

## 
##  Pearson's product-moment correlation
## 
## data:  table_travel_user$travel_time_mins and table_travel_user$user_type
## t = 44.664, df = 60437, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1710243 0.1864597
## sample estimates:
##      cor 
## 0.178753

#filtramos para los valores 1 y 2, eliminando trayectos de la empresa
z <- table_travel_user %>%
  filter(user_type=='1'|user_type=='2')

cor.test(z$travel_time_mins, z$user_type, method = c("pearson"))

## 
##  Pearson's product-moment correlation
## 
## data:  z$travel_time_mins and z$user_type
## t = 86.378, df = 54738, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.3389511 0.3536958
## sample estimates:
##       cor 
## 0.3463448

rm(table_travel_user, z)

No hay correlación entre el tipo de usuario y el tiempo. Se puede ver que eliminado los trayectos de los trabajadores, subre ligeramente el coeficiente de pearson, ya que los trayecto son más duraderos.

6.5 Análisis de correlación - user_age vs travel_time_mins

#Extraemos la variable distrito inicial eliminado Retiro 
table_travel_age <- data %>% 
  dplyr::select(c(travel_time_mins, user_age_range)) %>% 
  mutate(user_age_range = as.numeric(user_age_range))

#calculamos coeficiente de pearson
cor.test(table_travel_age$user_age_range, table_travel_age$travel_time_mins, method = c("pearson"))

## 
##  Pearson's product-moment correlation
## 
## data:  table_travel_age$user_age_range and table_travel_age$travel_time_mins
## t = -5.1455, df = 60437, p-value = 2.676e-07
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.02889337 -0.01295548
## sample estimates:
##         cor 
## -0.02092575

Al filtrar los usuarios de edad desconocidos el p_value pasa de negativo a ser >0,5, aunque el test no refleja correlación directa, si que parece que a mayor edad menor tiempo

#filtramos los desconocidos
table_trave_age_filter <- table_travel_age%>%
  filter(user_age_range!=0)

cor.test(table_trave_age_filter$user_age_range, table_trave_age_filter$travel_time_mins, method = c("pearson"))

## 
##  Pearson's product-moment correlation
## 
## data:  table_trave_age_filter$user_age_range and table_trave_age_filter$travel_time_mins
## t = 0.51208, df = 34963, p-value = 0.6086
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.007743348  0.013220007
## sample estimates:
##        cor 
## 0.00273863

rm(table_travel_age, table_trave_age_filter)

6.6 Análisis de correlación - hour.temp max vs unplug_hour

Estudio por hora de la posible correlación entre la hora de temperatura máxima e inicio del trayecto.

El test de pearson no de fiabilidad de correlación.
El dataset de Meteo con más meses si que podría haber sido de más utilidad.
También recogíamos la variable lluvia, pero en Agosto de 2018 la pluviometría fue de 0,6 L.

#Generamos la tabla con dos variables
table_hour_max <- data %>% 
  dplyr::select(c(unplug_hour, temp_max_hora)) %>% 
  mutate (unplug_hour= as.numeric(unplug_hour),
          temp_max_hora = as.numeric(temp_max_hora))


#Realizamos también el Test de Pearson
cor.test(table_hour_max$unplug_hour, table_hour_max$temp_max_hora, method = c("pearson"))

## 
##  Pearson's product-moment correlation
## 
## data:  table_hour_max$unplug_hour and table_hour_max$temp_max_hora
## t = -4.5098, df = 60437, p-value = 6.501e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.02630995 -0.01037044
## sample estimates:
##         cor 
## -0.01834136

rm(table_hour_max)

7 MODELOS DE CLUSTERING

7.1 PAM

Se utiliza PAM (Partición Alrededor de Medoids) que es un algoritmo de clasificación del tipo k-medoids. Escoge datapoints como centros y trabaja con una métrica arbitraria de distancias entre datapoints.

vars_cluster <- data %>% 
  filter(between (unplug_date,'01/08/2018','28/08/2018')) %>%
  filter( user_type == '1') %>%
  select(idunplug_station_ini,cp_ini,idplug_station_fin,cp_fin,travel_time,user_type,
         user_age_range,unplug_hour,distance,day, bank_holiday_flg) %>%
  mutate(
    day = case_when(
      day =='L' ~ 0,
      day =='M' ~ 1,
      day =='X' ~ 2,
      day =='J' ~ 3,
      day =='V' ~ 4,
      day =='S' ~ 5,
      day =='D' ~ 6),
    user_age_range=as.factor(user_age_range),
    unplug_hour=as.factor(unplug_hour),
    user_type=as.factor(user_type),
    idunplug_station_ini=as.factor(idunplug_station_ini),
    idplug_station_fin=as.factor(idplug_station_fin),
    cp_ini=as.factor(cp_ini),
    cp_fin=as.factor(cp_fin),
    bank_holiday_flg=as.factor(bank_holiday_flg))
   #filtramos para quedarnos con los tipos 1 y 2
   
#Selecionamos variables
data_cluster <- vars_cluster %>%
  select(idunplug_station_ini,cp_ini,idplug_station_fin,cp_fin,travel_time,user_type,
         user_age_range,unplug_hour,distance,day,bank_holiday_flg)

rm (vars_cluster)

7.1.1 CON OUTLIERS

En esta ejecucion del algoritmo se utilizan los datos con outliers (sin filtrado previo). Hacemos un subconjunto de 2500 muestras y calculamos el número de clusters recomendados para el Kmeans con el método del codo.

#Subsetting para clasterizar.
data_cluster_sample <- sample_n(data_cluster, 2500)

#Num cluster - Elbow method for KMEANS
fviz_nbclust(x = data_cluster_sample, FUNcluster = kmeans, method = "wss", k.max = 20) +
  geom_vline(xintercept = 5, linetype = 2)

Calculamos la distancia de Gower más representativa para el algoritmo PAM y utilizamos el método Silhouette, donde se escoge el primer pico cómo número posible de clusters.

data_cluster <- data_cluster %>%
  select(idunplug_station_ini,cp_ini,idplug_station_fin,cp_fin,travel_time,user_type,
         user_age_range,unplug_hour,distance,day)

rm (vars_cluster)

## Warning in rm(vars_cluster): object 'vars_cluster' not found

# PAM - Gower Distance
gower_dist <- daisy(data_cluster_sample, metric = "gower")
#2 clusters has the highest silhouette width. Let’s pick k = 5

#search for a number of clusters
sil_width <- c(NA)
for(i in 2:8){  
  pam_fit <- pam(gower_dist, diss = TRUE, k = i)  
  sil_width[i] <- pam_fit$silinfo$avg.width  
}
plot(1:8, sil_width,
     xlab = "Number of clusters",
     ylab = "Silhouette Width")
lines(1:8, sil_width)

Ejecutamos el algoritmo PAM para 4 clusteres

#executing the algorithm with 4 clusters
k <- 4
pam_fit <- pam(gower_dist, diss = TRUE, k)
pam_results <- data_cluster_sample %>%
  mutate(cluster = pam_fit$clustering) %>%
  group_by(cluster) %>%
  do(the_summary = summary(.))

#Visualization in a lower dimensional space
tsne_obj <- Rtsne(gower_dist, is_distance = TRUE)
tsne_data <- tsne_obj$Y %>%
  data.frame() %>%
  setNames(c("X", "Y")) %>%
  mutate(cluster = factor(pam_fit$clustering))
ggplot(aes(x = X, y = Y), data = tsne_data) +
  geom_point(aes(color = cluster))

Se observa que la clasterización no es muy buena con 4 clusters. Probamos con 5 que es lo que nos recomiendan ambos métodos de distancia.

Se obtiene una mejor representación, aunque debido a la muestra de 2500 muestras a veces la clasificación no es muy adecuada..

#The results looks like 3-4cluster, so we enter again.

k <- 5
pam_fit <- pam(gower_dist, diss = TRUE, k)
pam_results <- data_cluster_sample %>%
  mutate(cluster = pam_fit$clustering) %>%
  group_by(cluster) %>%
  do(the_summary = summary(.))

Visualización del modelo con 5 clusters.

#Visualization in a lower dimensional space
tsne_obj <- Rtsne(gower_dist, is_distance = TRUE)
tsne_data <- tsne_obj$Y %>%
  data.frame() %>%
  setNames(c("X", "Y")) %>%
  mutate(cluster = factor(pam_fit$clustering))
ggplot(aes(x = X, y = Y), data = tsne_data) +
  geom_point(aes(color = cluster))

Resumen de resultados los 5 clusters.

#Clusters’ scorcard
pam_results$the_summary

## [[1]]
##  idunplug_station_ini     cp_ini    idplug_station_fin     cp_fin   
##  80     : 30          28007  :102   64     : 96        28001  :108  
##  74     : 27          28009  : 62   91     : 17        28014  : 55  
##  86     : 19          28014  : 60   73     : 16        28009  : 54  
##  84     : 17          28012  : 55   100    : 15        28007  : 53  
##  73     : 16          28004  : 35   80     : 14        28012  : 49  
##  102    : 11          28006  : 24   74     : 13        28006  : 31  
##  (Other):317          (Other): 99   (Other):266        (Other): 87  
##   travel_time    user_type user_age_range  unplug_hour     distance   
##  Min.   :  157   1:437     0:281          18     : 46   Min.   :   0  
##  1st Qu.:  477             1:  0          19     : 35   1st Qu.:1231  
##  Median :  664             2:  3          20     : 34   Median :1894  
##  Mean   : 1038             3: 37          11     : 26   Mean   :1883  
##  3rd Qu.: 1026             4: 71          17     : 26   3rd Qu.:2391  
##  Max.   :20611             5: 43          22     : 24   Max.   :5661  
##                            6:  2          (Other):246                 
##       day       bank_holiday_flg    cluster 
##  Min.   :0.00   0: 77            Min.   :1  
##  1st Qu.:5.00   1:360            1st Qu.:1  
##  Median :5.00                    Median :1  
##  Mean   :4.76                    Mean   :1  
##  3rd Qu.:6.00                    3rd Qu.:1  
##  Max.   :6.00                    Max.   :1  
##                                             
## 
## [[2]]
##  idunplug_station_ini     cp_ini    idplug_station_fin     cp_fin   
##  90     :118          28009  :346   99     : 51        28004  :138  
##  65     : 46          28014  : 29   73     : 14        28012  : 55  
##  75     : 40          28004  : 28   86     : 14        28014  : 46  
##  66     : 35          28006  : 19   9      : 12        28006  : 42  
##  102    : 29          28012  : 13   170    : 11        28009  : 40  
##  101    : 27          28045  :  7   31     : 11        28007  : 31  
##  (Other):177          (Other): 30   (Other):359        (Other):120  
##   travel_time     user_type user_age_range  unplug_hour     distance   
##  Min.   : 116.0   1:472     0:119          15     : 81   Min.   :   0  
##  1st Qu.: 387.2             1:  1          9      : 45   1st Qu.:1081  
##  Median : 557.5             2:  1          8      : 34   Median :1694  
##  Mean   : 761.8             3: 27          14     : 31   Mean   :1748  
##  3rd Qu.: 835.8             4:242          20     : 30   3rd Qu.:2281  
##  Max.   :7568.0             5: 80          18     : 28   Max.   :5141  
##                             6:  2          (Other):223                 
##       day        bank_holiday_flg    cluster 
##  Min.   :0.000   0:427            Min.   :2  
##  1st Qu.:1.000   1: 45            1st Qu.:2  
##  Median :2.000                    Median :2  
##  Mean   :2.176                    Mean   :2  
##  3rd Qu.:3.000                    3rd Qu.:2  
##  Max.   :6.000                    Max.   :2  
##                                              
## 
## [[3]]
##  idunplug_station_ini     cp_ini    idplug_station_fin     cp_fin   
##  79     : 76          28007  :236   90     :111        28009  :309  
##  74     : 49          28012  : 67   65     : 39        28012  : 52  
##  80     : 34          28014  : 48   66     : 38        28007  : 48  
##  78     : 32          28004  : 45   75     : 36        28006  : 41  
##  73     : 19          28006  : 45   102    : 25        28014  : 37  
##  86     : 16          28045  : 21   107    : 17        28001  : 31  
##  (Other):380          (Other):144   (Other):340        (Other): 88  
##   travel_time      user_type user_age_range  unplug_hour     distance   
##  Min.   :    1.0   1:606     0:128          17     : 63   Min.   :   0  
##  1st Qu.:  494.5             1:  0          8      : 60   1st Qu.:1502  
##  Median :  667.5             2:  3          7      : 48   Median :2013  
##  Mean   :  868.6             3: 41          19     : 40   Mean   :2090  
##  3rd Qu.:  893.5             4:339          20     : 38   3rd Qu.:2533  
##  Max.   :18543.0             5: 93          21     : 37   Max.   :5779  
##                              6:  2          (Other):320                 
##       day        bank_holiday_flg    cluster 
##  Min.   :0.000   0:564            Min.   :3  
##  1st Qu.:1.000   1: 42            1st Qu.:3  
##  Median :2.000                    Median :3  
##  Mean   :2.076                    Mean   :3  
##  3rd Qu.:3.000                    3rd Qu.:3  
##  Max.   :6.000                    Max.   :3  
##                                              
## 
## [[4]]
##  idunplug_station_ini     cp_ini    idplug_station_fin     cp_fin   
##  64     :147          28001  :172   84     : 80        28014  :236  
##  73     : 19          28012  : 60   91     : 43        28009  : 59  
##  99     : 19          28009  : 59   85     : 28        28007  : 47  
##  65     : 15          28006  : 42   73     : 27        28006  : 39  
##  74     : 15          28014  : 39   69     : 21        28012  : 31  
##  86     : 13          28004  : 36   72     : 18        28016  : 28  
##  (Other):326          (Other):146   (Other):337        (Other):114  
##   travel_time      user_type user_age_range  unplug_hour     distance   
##  Min.   :  113.0   1:554     0:369          13     : 51   Min.   :   0  
##  1st Qu.:  428.0             1:  0          19     : 45   1st Qu.:1172  
##  Median :  610.0             2:  3          20     : 39   Median :1693  
##  Mean   :  891.0             3: 26          15     : 35   Mean   :1770  
##  3rd Qu.:  907.8             4: 92          7      : 35   3rd Qu.:2217  
##  Max.   :17249.0             5: 60          18     : 32   Max.   :5618  
##                              6:  4          (Other):317                 
##       day        bank_holiday_flg    cluster 
##  Min.   :0.000   0:497            Min.   :4  
##  1st Qu.:2.000   1: 57            1st Qu.:4  
##  Median :3.000                    Median :4  
##  Mean   :2.684                    Mean   :4  
##  3rd Qu.:4.000                    3rd Qu.:4  
##  Max.   :6.000                    Max.   :4  
##                                              
## 
## [[5]]
##  idunplug_station_ini     cp_ini    idplug_station_fin     cp_fin   
##  91     : 57          28014  :162   74     : 86        28007  :212  
##  84     : 32          28012  : 42   79     : 32        28009  : 36  
##  69     : 26          28009  : 41   78     : 28        28012  : 33  
##  99     : 18          28004  : 38   80     : 24        28006  : 31  
##  85     : 16          28007  : 32   100    : 11        28014  : 18  
##  74     : 12          28006  : 30   73     : 11        28001  : 16  
##  (Other):270          (Other): 86   (Other):239        (Other): 85  
##   travel_time      user_type user_age_range  unplug_hour     distance   
##  Min.   :  115.0   1:431     0:105          14     : 75   Min.   :   0  
##  1st Qu.:  466.5             1:  3          15     : 31   1st Qu.:1299  
##  Median :  654.0             2:  4          19     : 29   Median :1862  
##  Mean   :  872.2             3: 18          20     : 29   Mean   :1976  
##  3rd Qu.:  948.0             4: 64          8      : 27   3rd Qu.:2420  
##  Max.   :10479.0             5:233          16     : 26   Max.   :5552  
##                              6:  4          (Other):214                 
##       day        bank_holiday_flg    cluster 
##  Min.   :0.000   0:386            Min.   :5  
##  1st Qu.:0.000   1: 45            1st Qu.:5  
##  Median :1.000                    Median :5  
##  Mean   :1.838                    Mean   :5  
##  3rd Qu.:3.000                    3rd Qu.:5  
##  Max.   :6.000                    Max.   :5  
##

rm(pam_results)
rm(data_cluster_sample)

La parte positiva del modelo de PAM es que obtenemos unos resultados que aunque habría que ejecutar con una muestra más amplia o mejor definida, nos da unos patrones que nos permitirían explicar de manera clara a negocio los diferentes clusters.

7.2 CLARA

Se utiliza CLARA (Clustering Large Applications) es un algoritmo de clasificacion basado en K-Means que utiliza tecnicaas de muestreo para caracterizar grandes datasets. Toma un numero pre-determinado de clusters.

# Seleccion de variables para el modelo
vars_modelo <- c('idunplug_station_ini','cp_ini','idplug_station_fin','cp_fin','travel_time','user_type', 'user_age_range','unplug_hour','day','distance')

7.2.1 CON OUTLIERS

En esta ejecucion del algoritmo se utilizan los datos con outliers (sin filtrado previo) y 6 clusters.

# CON OUTLIERS

# dataset
data_modelo <- data %>% 
  select(vars_modelo)

# Convertir data types
data_modelo[sapply(data_modelo, is.character)] <- lapply(data_modelo[sapply(data_modelo, is.character)], as.factor)


# Compute CLARA
# clara_result <- clara(data_modelo, 10, samples = 2500, pamLike = TRUE)
load(file = 'clara_result.RData')

# Add the point classifications to the original data, use this:
data_modelo_clusters <- cbind(data_modelo, cluster = clara_result$cluster)
head(data_modelo_clusters, n = 5)

##   idunplug_station_ini cp_ini idplug_station_fin cp_fin travel_time
## 1                   99  28004                103  28009         291
## 2                   86  28012                 77  28009         549
## 3                  100  28006                 52  28012         632
## 4                   44  28012                102  28009         847
## 5                  165  28046                 91  28014        1289
##   user_type user_age_range unplug_hour day distance cluster
## 1         1              5           1   J 1334.102       1
## 2         1              4           1   J 1970.717       2
## 3         1              0           1   J 2293.919       3
## 4         1              0           1   J 2319.953       3
## 5         1              0           1   J 1902.985       3

rm(data_modelo_clusters)

En el grafico se observa que el cluster-1 agrupa todos los aoutliers del dataset.

# Cluster plot
fviz_cluster(clara_result, stand = FALSE, geom = "point", pointsize = 1)

7.2.2 SIN OUTLIERS

En esta ejecucion del algortimo se toma el datatset original filtrado por quartiles extremos con la funcion boxplot.stats()$out. Se toman 5 clusters.

# SIN OUTLIERS

# datatset
data_modelo_nooutliers <- data_no_outliers %>% 
  select(vars_modelo)

# Convertir data types
data_modelo_nooutliers[sapply(data_modelo_nooutliers, is.character)] <- lapply(data_modelo_nooutliers[sapply(data_modelo_nooutliers, is.character)], as.factor)

# Compute CLARA
# clara_result_nooutliers <- clara(data_modelo_nooutliers, 10, samples = 2500, pamLike = TRUE)
load(file = 'clara_result_nooutliers.RData')

# Add the point classifications to the original data, use this:
data_modelo_nooutliers_clusters <- cbind(data_modelo_nooutliers, cluster = clara_result_nooutliers$cluster)
head(data_modelo_nooutliers_clusters, n = 4)

##   idunplug_station_ini cp_ini idplug_station_fin cp_fin travel_time
## 1                   99  28004                103  28009         291
## 2                   86  28012                 77  28009         549
## 3                  100  28006                 52  28012         632
## 4                   44  28012                102  28009         847
##   user_type user_age_range unplug_hour day distance cluster
## 1         1              5           1   J 1334.102       1
## 2         1              4           1   J 1970.717       2
## 3         1              0           1   J 2293.919       2
## 4         1              0           1   J 2319.953       3

rm(data_modelo_nooutliers_clusters)

Se observa que los grupos est’an agrupados entre si, pero no hay mucha distancia entre clusters.

# Cluster plot
fviz_cluster(clara_result_nooutliers, stand = FALSE, geom = "point", pointsize = 1)

8 CONCLUSIONES

-Los picos de uso de Bicimad se observan muy marcados con el inicio y fin de la jornada laboral, podría ser interesantes realizar acciones comerciales dirigidas hacia esas horas

-Excluyendo los de edad desconocida, el 86% de usuarios son adultos (27-65 años)

-Uso constante casi rutinario a nivel de semana, cayendo en los fines de semana, más pronunciado esto en los usuarios de edad adulta (27-65 años)

-En lo referente a los clusters, aunque no se observan grandes diferencias entre los centroides, si aparecen bien definidos y explicados por las variables seleccionadas

Comparando con los datos de Enero: - El uso del servicio baja en las horas de mayor calor 15-18 horas, se intensifica entre las 7-9 de la mañana. - En verano el uso entre las 19-22 horas se amplía - Los trayectos son más cortos, destacan muchos movimiento en el propio retiro o con estaciones centrales.

PROXIMOS PASOS

-Utilizar herramientas desktop para no tener problemáticas de memoria

-Cruzar con datasets de eventos ligado a zonas

-Volver a ejecutar el proyecto con datos de otros periodos, pues están algo sesgados por la temporalidad, para poder ver como se adaptan.

Modulo5_grupal - Análisis de los datos de Bicimad mediante clustering

Alejandra Mendez, Alvaro Cerrato, Francisco Martin, Juan Antonio Castro

9/6/2019

1 SYSTEM SETUP

1.1 Paquetes de R

2 DATA WRANGLING & FEATURE ENGINEERING

2.1 Tablas auxiliares

3 ANALISIS EXPLORATORIO

3.0.1 Distrito de inicio

3.0.2 Codigo postal de inicio

3.0.3 Estacion de inicio

3.0.4 Distrito de fin

3.0.5 Codigo postal de fin

3.0.6 Estacion de fin

3.0.7 Tipo de usuario

3.0.8 Rango de edad de usuario

3.0.9 Codigo postal de usuario

3.0.10 Dia de la semana del trayecto

3.0.11 Festivos/Laborables

3.0.12 Hora de inicio del trayecto

3.1 ESTADISTICOS PRINCIPALES

3.1.1 Duracion del trayecto

3.1.2 Distancia del trayecto

3.1.3 Temperatura maxima

3.1.4 Temperatura minima

3.1.5 Velocidad del viento maximas

4 ANALISIS DE OUTLIERS

4.1 TRAVEL_TIME

4.2 DISTANCE

5 ANALISIS BIVARIANTE

5.1 Tb.Contingencia - Edad de usuario v. Tipo de usuario

5.2 Estudio por Edad de los trayectos con distrito final enel Retiro y distrito inicial diferente

5.3 Estudio por Edad de los trayectos con origen distrito Retiro y distrito final diferente al Retiro.

6 ANALISIS DE CORRELACION

6.1 Matriz de correlación - Pearson

6.2 Matriz de correlación - Spearman

6.3 Análisis gráfico - distance vs travel_time_mins

6.4 Análisis de correlación - user_type vs travel_time_mins

6.5 Análisis de correlación - user_age vs travel_time_mins

6.6 Análisis de correlación - hour.temp max vs unplug_hour

7 MODELOS DE CLUSTERING

7.1 PAM

7.1.1 CON OUTLIERS

7.2 CLARA

7.2.1 CON OUTLIERS

7.2.2 SIN OUTLIERS

8 CONCLUSIONES